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ABSTRACT 

We investigate the chemistry of ion molecules in protoplanetary disks, mo¬ 
tivated by the detection of N2H+ ring around TW Hya. While the ring inner 
radius coincides with the CO snow line, it is not apparent why N2H''' is abun¬ 
dant outside the CO snow line in spite of the similar sublimation temperatures 
of CO and N2. Using the full gas-grain network model, we reproduced the N2H'*' 
ring in a disk model with millimeter grains. The chemical conversion of CO and 
N2 to less volatile species (sink effect hereinafter) is found to affect the N2H''' 
distribution. Since the efficiency of the sink depends on various parameters such 
as activation barriers of grain surface reactions, which are not well constrained, 
we also constructed the no-sink model; the total (gas and ice) CO and N2 abun¬ 
dances are set constant, and their gaseous abundances are given by the balance 
between adsorption and desorption. Abundances of molecular ions in the no-sink 
model are calculated by analytical formulas, which are derived by analyzing the 
full-network model. The N2H’'' ring is reproduced by the no-sink model, as well. 
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The 2D (R-Z) distribution of N 2 H+, however, is different among the full-network 
model and no-sink model. The column density of N 2 H+ in the no-sink model 
depends sensitively on the desorption rate of CO and N 2 , and the flux of cosmic 
ray. We also found that N 2 H+ abundance can peak at the temperature slightly 
below the CO sublimation, even if the desorption energies of CO and N 2 are the 
same. 

Subject headings: astrochemistry — protoplanetary disks 


Introduction 


A ring of N 2 H’'‘ emission was found in th e disk around T W Hya using ALMA (Atacama 
Large Millimeter and sub millimeter Array) fjOi et al.l 120131) . N 2 H+ is considered to be a 
useful probe of CO snow line, because it is destroyed by the proton transfer to CO 


N 2 H+ + CO ^ HCO^ 


Ns 


( 1 ) 


The inner radius of N 2 H'^ ring is indeed consistent with the CO snow line predicted from a 
disk model of TW Hya. The anti-correlatio n of N 2 H’'' and CO emission is also often observed 
in prestellar cores (e.g. iTafalla et al.l 12004 1. which strengthen the above statement. 


The bright emission of N 2 H+ outside the CO snow line is, however, puzzli ng. Laboratory 


experiments show that t he sublimation temperatures of CO and N 2 are similar (jCollings et ah 


2004 lOberg et al.ll2005l) . Since N 2 H''' is formed by protonation of N 2 , N 2 H''' should also be 
depleted outside the CO snow line. In the case of prestellar cores, there are two possible 
explanations for the survival of N 2 H’*’ in the CO depleted region. Firstly, destruction of 
N 2 H''‘ is temporally suppressed by the freeze-out of its major reactant, CO. Secondly, the 
How formation of No from N atom in molecular clouds helps the temporal survival of N 2 
f Aikawa et al. 2001 : Maret et ah 2006 : Bergin &: Tafalla 2007 ). The hrst mechanism could 
work in the disk as well, in a narrow range of temperatures where the positive effect of 
depletion on N 2 H+ abundance (i.e. decrease of CO) wins against the negative effect (i.e. 
depletion of N 2 ). The second mechanism, on the other hand, seems irrelevant for disks. In 
protoplanetary disks, in which the gas density is much higher than molecular clouds, it is 
not likely that N atom is more abundant than N 2 except for the photodissociation layer at 
disk surface. 

In chemistry models of disks, it is often found that CO is depleted even in regions 
warmer than its sublimation ternperature 20 K) via conve rsion to less volatile species 
flAikawa et al.lll997l: iBergin et al.ll2014J: iFuruva fc Aikawal 120141) . Since the conversion works 
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as a sink in th e chemical r e action network of gaseous species, we call it the sink effect in the 
present work. iFavre et ahl (120131) showed, for the first tirn e, observational evidences for such 
CO depletion towards TW Hya. iFuruva fc Aikawal (120141) showed that N 2 is also subject to 
the sink effect; it is converted to NH3 via the gas-phase reactions of 


N 2 + H+ ^ N 2 H+ + H 2 
N 2 H+ + e ^ NH + N 


( 2 ) 

( 3 ) 


and subsequent hydrogenations of NH on grain surfaces. It should be noted that this conver¬ 
sion is not efficient when CO is abundant, because N 2 H+ mainly reacts with CO to reform N 2 . 
The longer timescale required for N 2 depletion could explain the N 2 H+ ring. But then the 
lifetime of the ring should be similar to the difference in the depletion ti mescales of CO and 
N 2 , w hich is rather short < 10® yr, if the ionization rate is ~ 5 x 10“^^ s“^ (IFuruva fc Aikawa 
20141 ). Since N 2 is considered to be a major nitrogen carrier in protoplanetary disks, but 
cannot be directly observed, it is important to constrain its distribution based on N 2 H''' ob¬ 
servation. The formation mechanism of N 2 H''' ring is thus worth investigation via theoretical 
calculations. 


The abundance of N 2 H''' could also be a probe of ionization degree in the regions of CO 
depletion. The ionization degree depends on the gas density and ionization rate, which is an 
important parameter for both chemical and physical evolutions of the disk. Since ionizations 
trigger the chemical reactions, as in molecular clou ds, the timescale of chemical evolution 
depends on the ionization rate (lAikawa et al.lll997l) . The ionization degree determines the 
coupli ng with magnetic felds, e.g. which disk region is subject to magneto-rotational insta¬ 
bility (iBalbus fc Hawlevlll99l[) . The major io nization sources ye X-ray from the central star. 


cosni i c ray, and decay of radioactive nuclei (jClassgold et ahl Il997l: lUmebavashi fc Nakano 


19881: ICleeves et ah 1120131) . The ionization rate thus depends on the flux and hardness of 


X-ray radiation and the abundances of radioactive nuclei. In addition, the stellar winds 
and/or magnetic helds of the star-disk system could prevent the penetration of the cosmic 
ray to the disk. Since these parameters are unknown and could vary among objects, the 
ionization degree should be probed via observations. The observation of ionization degree is 
principally possible using major ion molecules. But a quantitative estimate of the ionization 
degree from the observational data is not straightforward, because the g as density, maior ion 


molecules^ and their abundances ch ange spatially within the disk (e.g lAikawa et ah 


2002 


Qi et ah II 2 OO 8 I: ICleeves et ah 1120141) . 


In this work, we investigate the spatial distribution of N 2 H+ and other major molecular 
ions in protoplanetary disks. Firstly, we calculate the full chemical network model, which 
includes gas-phase and grain-surface reactions. The radial distribution of N 2 H’^ column 
density has a peak around the radius of CO sublimation temperature; i.e. the N 2 H+ ring 




































4 


is reproduced. We found that the sink effect on CO signihcantly affects the distribution of 
N 2 H+. By analyzing the full chemical network, we also found that the abundances of major 
molecular ions, HCO^ and N 2 H’'' can be described by analytical functions of gas density, 
temperature, ionization rate, and abundances of CO and N 2 . Such formulas are useful in 
deriving the ionization degree from the observations of molecular ions, combined with the 
dust continuum and CO observations to constrain the gas density and CO abundance. 

While the sink effect plays an important role in determining the N 2 H"*^ distribution in 
the full network model, the efficiency of the sink depends on various parameters which are 
not well constrained yet. For example, the conversion timescale of CO and N 2 to less volatile 
species depends on ionization rate and activation barriers of the reactions. The sink effect 
could be less significant, if the disk is turbulent and the diffusion timescale is shorter than the 
conversion timescale. Since there are various chemical paths for the conversion, it is difficult 
to directly control the efficiency of sink in the full-network model. In order to investigate the 
abundances of N 2 H"'' and other molecular ions in the limit of no sink, we construct ”no-sink” 
model; we assume that the total (gas and ice) abundances of CO and N 2 are constant, and 
that their gas-phase abundances are given by the balance between adsorption and desorption. 
The abundances of molecular ions are calculated using the analytical formulas. 

The plan of this paper is as follows. In §2, we describe our disk model and chemical 
model. The molecular distributions in the full-network model are presented in section §3. 
In §4, we derive analytical formulas for abundances of electron, H^, HCO’*' and N 2 II+. The 
formulas are used to calculate their abundances in the disk model, which are compared with 
the results of full-network calculation. Section 5 presents the results of no-sink model. Since 
the no-sink model is analytical, we can easily investigate the dependence of N 2 H''’ column 
density on desorption rate (i.e. sublimation temperature) of CO and N 2 . We also investigate 
how much the N 2 II’'' column density is reduced, if the penetration of cosmic rays to the 
disk is hampered by the stellar wind and/or magnetic fields. We summarize our results and 
conclusions in 56. 


2. Models 


We adopt the same disk model as in iFuruva fc Aikawal (1201^. In addition to the disk 
model with small grains, which is assumed in lFuru~ fc Aikawal J2OI4I) . we also investigate 
a disk model in which the dust grains have grown up to a radius of 1 mm. The turbulent 
mixing is not explic itly included in the prese nt work. The model of full chemical network is 
also the same as in iFuruva fc Aikawal (120141) . but includes some updates. The models are 
briefly described in the following. 
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2.1. Disk Model 


Since we aim to understand the mechanism to form N 2 H”^ ring, rather than to investi¬ 
gate the disk structure of a specific object (i.e. TW Hya), we adopt a steady, axisymmetric 
Keplerian disk around a T Tauri star. The stellar mass, radius and effective temperature 
are M* = O.SMq, = 2Rq and T* = 4000 K, respectively. The disk structure is given 
by solving the radiation transfer, thermal balance of gas and dust, and hydrostatic equi¬ 
librium in the v ertical direction in th e disk. Basic equations and calculation procedures 
are described in iNomura et ah (120071) . We assume the stellar UV and X-ray luminosity 
of 10^^ erg s“^ and ip^o res pectively. The cosmic ray ionization rate of H 2 is set 

to be 5 X 10“^^ s“^ ( Dalgarno 20061). while the ionization rat e by the decay of radioactive 
nuclei is set to be 1 x 10“^® s“^ fjUmebavashi fc Nakano II2009I) . The dust-to-gas mass ratio 
is 0.01. We consider two disk models: one with the dark cloud dust and the other with 


m illimeter-sized grains (lAikawa &: Nomura I l2006h . The former assumes the dust properties 
of IWeingartner &: Draine (120011) [R^ = 5.5, be = 3 x 10“^, case B); while the silicate grains 
have a rather steep size distribution with the maximum radius of ~ 0.2 /rm, the carbona¬ 
ceous grains have PAH-like properties in the small-size limit and graphite-like properties at 
larger sizes. The maximum size of carbonaceous grain is ~ 10 /im. In the latter model, 
we assume the power-law size distribution of d ust grains dn(a)/da (x where a is the 

grain radius, referring to the ISM dust model of iMathis et ah (jl977f ). but the minimum and 
maximum sizes are set to be 0.01 /im and 1 mm, respectively. It would be more appropriate 
for T Tauri disks than the dark cloud dust model, since the grain growth is indicated by 
the disk observations (e.g. IWilliams &: Cieza 1120111) . The dust opacities for the two mod¬ 
els are calculated using Mie theory. The gas temperature, dust temperature, and density 
distributions in the disk are calculated self-consistently, by considering various heating and 
cooling mechanisms. Figure [T] shows the distribution of gas density, gas temperature, dust 
temperature, and ionization rate by X-ray in our models. Cosmic ray ionization dominates 
in the midplane where the X-ray ionization rate is < 5 x 10“^^ s“^. The temperature in 
the model with millimeter grains is lower than that in the dark cloud dust model due to 
the lower dust opacity at a given disk height {Z) to receive the radiation from the central 
star. Column density in these disk models are determined by assuming a steady state disk 
structure with constant viscosity and accretion rate (although we consider turbulent diffu¬ 
sion and/or radial accretion only implicitly in the chemical model), so that the masses of the 
two disks are slightly different: 1 x 10 “^Mq for the dark cloud dust model and 1.7 x 10“^Mp 
for the millimeter grain model. 
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2.2. Chemical Model: Full Network 


Our chemical network is based on iGarrod fc Herbst (120061) . We added photo-ionization, 
photo-dissociation and photo-desorption by UV radiat ion from the central star, self-shielding 


of H 2 , CO and N 2 flFuruva et ahll2013l: iLi et ah 1120131) . X-ray chemistry and charge balance 


of dust grains. Although our model includes Deuteration and ortho/para states of several 
species such as H 2 and H 3 (Hincelin et ah in prep, Furuya et ah in prep), we present only the 
molecular abundances, i.e. the sum of isotopomers and o/p states. The D/H ratio and o/p 
ratios will be presented in forthcoming papers. Our model consists of two phases, gas phase 
and ice mantle; i.e. we do not discriminate layers of ice mantles, unl ess otherwise stated. 


Desor ption energies (Fjes) of assorted species are listed in Table 1 in iFuruva &: Aikawa 


(1201411 . Desorption energies of atomic hydrogen, CO and N 2 are s et to be 600 K, 1150 K, and 


1000 K, respectively; they are the values on water ice substrates (lAl-Halabi fc van Dishoeck 


2 OO 7 I: ICarrod &: Herbst II2006I) . The sublimation temperature of CO and N 2 are then ~ 23 
K and ~ 19 K, when the gas density is 10® cm“^. We investigate the dependence of N 2 H’'' 
abundance on the desorption energies of CO and N 2 in §3 and §5. 

Adsorbed species on grains migrate via , thermal hopping and rea ct with each other when 
they meet. We adopt the modihed rate of ICaselli et al.l (jl998L 1200211 for grain surface reac¬ 
tions of H atom. The adsorption rate of gaseous species onto grain surfaces and grain surface 
reaction rates (e.g. if the rate is limited by the accretion of gaseous particle) depend on the 
size distribution of gra ins. Ideally, the grain size distribution should be taken into account 
(lAcharvva et al.l 1201 ll ). Most chemical models of disks and molecular clouds, however, as¬ 
sume a single size of 0.1 pm, for simplicity, which we follow in the present work. The rate 
coefficients of gas-dust interactions (e.g. adsorption) are basically proportional to the total 
grain surface area. The assumption of a single grain size in the chemical model is thus a 
reasonable approximation, as long as the total surface area of grains are consistent with that 
in the physical disk model. The total surface area of the 0.1 pm dust model agrees with that 
of our dark cloud dust model within a factor of a few. In the chemical model for the disk 
model with millimeter grains, we adopt the same uniform grain size (0.1 pm), but decrease 
the dust-gas ratio by one order of magnitude; according to the power-law size distribution, 
the number of small grains, which dominate in the grain sur face area, is decreased corn pared 
with the dark cloud dust model by an order of magnitude. (lAikawa &: Nomura 1120061) . One 
caveat for this single-size approximation is that we may underestimate the rate of grain 
surface recombination. In the dense regions, such as disk midplane at small radii, recom¬ 
binations are more efficient on grain surfaces than in the gas phase. Due to the Coulomb 
focusing, the cross section of grain-surface recombination is much larger than the geometrical 
cross section of dust grains. The rate of grain surface recombination thus cannot be scaled 
by the total grain surface area. 
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The photodissociation rates are calculated by convolving the attenuated stellar and 
interstellar UV spect rum and wavelength - dependent photodissociati on cross sections at each 


position in the disk flAikawa et ah 


2002 


van Dishoeck et ah 


2006h . UV radiation induced 

by X-ray and cosmic ray fiGredel et ah Il989l) is also taken into account. For ice mantle 
species, we assume that only the uppermost layers can be dissociated; i.e. while the UV 
radiation can penetrate into deep layers of ice mantle, we assume that the photo products in 
deeper layers recombine immediately. The effective rate of photodissociation is thus reduced. 
Considering the fluffiness and pores on the grain surfaces, we let the uppermost two layers, 
rather than one, to be dissociated. 

We take into account three non-thermal desorption processes: photodesorption, stochas¬ 
tic heating by cosmic-rays, and reactive desorption. We adopt the photodesorption yields 
per inci dent FUV photon d er ived from the labo r atory experiment fo r H 2 O, CO 2 , CO, O 2 , 
and N 2 f Oberg et ah 2009a b; Favolle et ah 2011 : Favolle et ah 2013). A yield is set to 10“^ 
for other species. 

Initial molecular abundance in the disk is given by ca l culati ng the mo lecular evolution 


in a s tar-forming core m odel of iMasunaga and Inutsuka (1200011 (see also iMasunaga et ah 


19981: lAikawa et al.l 1200811 . The initial abundances of assorted molecules and the elemental 
abundances in our model are listed in Table |2j Major carriers of oxygen and carbon are 
H 2 O and CO, while the major N-bearing species are NH 3 and N 2 . We adopt the low- 
metal abundance; i.e. the abundances of metals such as Mg and Si are about two orders of 
magnitude lower than observed in diffuse clouds. 

We calculate the chemical reaction network (i.e. rate equations) as an initial value 
problem at each position in the disk. As we will see in §3, the abundances of electron, 
HCO"*", N 2 H’'‘, and Hj]' reach the steady state, which are determined by the ionization rate, 
gas density, temperature, and the abundances of CO and N 2 , in a short timescale. On the 
other hand, CO and N 2 decrease slowly with time mainly due to the sink effect. Vertical 
diffusion and radial accretion, which are not explicitly i ncluded in the present work, could 
suppress or slow-down the sink effect flFuruva et al.ll2013l:lFuruva fc Aikawall2014fl . Therefore 
we present molecular abundances at an early time 1 x 10 ^ yr, as well as at the typical timescale 
of T Tauri stars rsj 10 ® yr. 


2.3. Chemical Model: No-Sink Model 

In the full-network model, we will see that the distribution of N 2 II'*’ is signihcantly af¬ 
fected by the CO depletion via sink effect. Efficiency of the sink effect, however, depends on 














































various parameters such as CO2 formation (C O + OH) rate on dust grains, initial CO ab un- 
dance, ionization rate, and turbulent mixing flBergin et al.ll2014j: iFuruva fc Aikawall2014l) . A 
strong vertical turbulence, for example, tends to smooth out the molecular abundances, so 
that the local CO abundance minima due to the sink could be less si gnificant. Althoug h the 
observation of TW Hya indicates the CO depletion via sink effect (IFavre et al.l 120131 1. the 
spatial distribution of CO abundance is not well constrained yet. It is therefore useful to 
calculate the distribution of N2H’'' in a model without the sink effect. 

In the no-sink model, we assume that the sum of gas-phase and ice-mantle abundances 
of CO is equal to its canonical abundance: i.e. 1 x 10“^ relative to the hydrogen nuclei. The 
gas-phase abundance of CO is given by a simple balance between adsorption and desorption; 


^COgas 

^COice 


z/exp(-%^) + z/rcRC'Feexp(-^%^) 


SttOj 


^^COgas + ^^COice = 10 


( 4 ) 

( 5 ) 


where ricogas and ricoice are number densities of CO in the gas phase and ice mantle. While 
the first term in the numerator represents the thermal desorption, the second term represents 
the non-thermal desorption. Although our full network model includes various mechanisms 
of non-thermal desorption, here we consider only the stochastic heating by cosmic-rays for 
simplicity. It is the effective deso rption mechanism in the cold midplane for species with 
relatively low desorption energies flHasegawa fc Herbstlll993l) . The frequency of CO oscilla¬ 
tion on grain surface u is set to be 10^^ s“^. When a cosmic-ray particle hits a dust grain, 
the grain is heated temporally for tcr ~ 10“^ sec. The peak temperature of the temporal 
heating is set to be T^ax = 70 K. The rate for a grain to encounter Fe ion particle, which is 
the most efficient in dust heating among cosmic-ray particles, is Cpe = 3 x 10“^^ s“^. The 
denominator represents the sticking rate of gaseous CO onto grain surfaces. The sticking 
probability on collision is set to be S' = 1.0. The grain size a is 0.1 /rm (see §2.2), and Uth is 
the thermal velocity of CO particle. The abundance of gaseous N2 is formulated similarly, 
with the total N2 abundance set to be 4.5 x 10“®. It should be noted that this desorption rate 
by the cosmic-ray heating is a rough estimate. In reality and in our physical disk models, the 
grains actually have a size d i stribu tion, and the parameters such as tcr and T^ax depend on 
the grain size. iLeger et ah fjl985l ). however, showed that the desorption rate via cosmic-ray 
heating does not sensitively depend on the grain size as long as the grains are small < 0.2 
fim. Since the small grains contribute most to the total surface area of grains, the desorption 
rate obtained here would be reasonable. 


In section §4, we analyze the results of full network model to find out that the abundances 
of electron, Hg , HCO'*', and N2H+ can be well described by analytical formulas, which are 
the functions of density, temperature, ionization rate, and abundances of CO and N2. We 
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use these analytical formulas to obtain the molecular ion abundances in the no-sink model. 
A combination of the analytical formulas of molecular ions and equilibrium abundances of 
gaseous CO and N2 (eq. S] and ED makes it very easy to investigate the dependence of N2H+ 
abundance on various parameters such as desorption energies of CO and N2 and ionization 
rate in the disk. 


3. Results: Full Network 
3.1. Disk with dark cloud dust 

Figure [2] shows the distributions of CO, HCO+, N2, N2H+, H3 and electron abundances 
at the time of 1 x 10^ and 9.3 x 10^ yr in the model with dark cloud dust. N2H'^ exists 
mostly in the upper layers of the disk. is more abundant than HCO’*' in the disk surface 
due to a relatively high abundance ratio of electron to CO (see §4.4). N2H+ is thus kept 
abundant there, via N2 -|- H^, in spite of the destruction via the reaction ([ID. 

The dashed lines in the panels of CO and N2 depict the position where the gas-phase 
abundance and ice abundance become equal in the adsorption-desorption equilibrium: 

^ ^ ^exp(-f^) + zyrcRCFeexp(-^) ^ ^ 

^ice StTO, Udust^th 

where Ugas and nice are number densities of CO (or N2) in the gas phase and ice mantle. It 
defines the ’’snow surface” of volatile species in the disk. We define the CO and N2 ’’snow 
line” as the radius at which the equation (ED is satisfied in the midplane. 

In Figure 121 we can see that both CO and N2 are depleted in the midplane even inside 
their snow lines, especially at the later stage, t = 9.3 x 10^ yr; CO is converted to CO2 
ice, CH3OH ice and hydrocarbons such as CH4 ice and C2H6 ice, while N2 is converted to 
NH3 ice. N2H+ is not abundant in the CO depleted region, since its mother molecule, N2, 
is depleted there as well. At f = 1 x 10^ yr, the sink effect is still moderate, and N2H+ 
abundance in the midplane has a local peak at ~ 200 AU. N2 is depleted in the outer radius, 
while N2H+ is destroyed by CO in the inner radius. 


3.2. Disk with millimeter grains 

Figure EJshows the distributions of CO, HCO+, N2, N2H+, H3 and electron abundances 
in the disk model with millimeter grains. The midplane temperature is warmer than 20 
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K and 30 K inside ~ 40 AU and ~ 10 AU, respectively. Althongh the sink effect is less 
significant dne to the smaller total snrface area of dust grains than in the model with dark 
cloud dust, CO is converted to CO 2 ice via the grain-surface reaction of CO -|- OH, and is 
depleted from the gas phase even in the intermediate {Z/R ~ 0.2) layers above the dashed 
line, where the dust temperature is higher than the sublimation temperature of CO. In the 
midplane region with T < 20 K, on the other hand, major carbon reservoirs are CO ice, CO 2 
ice and CH 3 OH ice. 

While the spatial distributions of CO and N 2 are similar in the dark cloud dust model, 
they are significantly different in the millimeter grain model; N 2 is abundant in a layer at 
ZjR ^ 0.2, where CO is depleted. In regions closer to the midplane, N 2 ice and NH 3 ice 
are the dominant N-bearing species. The conversion of N 2 to NH 3 ice proceeds via the 
gas-phase reactions of ([2]) and Q, and subsequent hydrogenation of NH on grain surfaces. 
At Z/i? ~ 0.2, the photodissociation of NH (NH —)■ N -|- H) is effective, which prevent the 
conversion of N 2 to NH 3 ice. The product of photodissociation, N atom, is converted back 
to N 2 via the reaction of N -|- NO. In other words, a deeper penetration of UV radiation in 
the millimeter grain model is a key to suppress the sink effect on N 2 . 

Another key difference between the models with dark cloud dust and millimeter grains is 
that Hj]" is not the dominant ion in the disk surface in the latter model. Since the UV shielding 
via dust grains is less effective in the millimeter grain model, UV radiation ionizes atoms, 
such as C and S, to make electrons abundant. Dissociative recombination of a molecular ion 
with an electron is much more efficient than the radiative recombination of an atomic ion. 
Thus the photoionization of atoms results in depletion of H^ and N 2 H’'' in the disk surface 
(e.g. Z/i? > 0.3). In lower layers, the distribution of N 2 H’'' is basically similar to that of N 2 , 
except for the midplane in inner radius where CO is abundant (i.e. R < a. few tens of AU 
at t = 1 X 10^ AU). 

It should be noted that the layer with abundant N 2 and CO depletion (and thus abun¬ 
dant N 2 H+) at ZjR ^ 0.2 is formed via the combination of sink effect on CO and photodis¬ 
sociation of NH, rather than via the difference in sublimation temperature of CO and N 2 . 
We performed a calculation of the same disk model but setting the both desorption energies 
of CO and N 2 to be 1000 K. Distributions of CO, HCO'*', N 2 and N 2 H+ are shown in Figure 
m The resultant distribution of molecules are basically the same as in Figure [3l 
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3.3. Column densities 


Although we do not aim to construct a best-fit model for a specific object, it is useful to 
briefly compare our models with observations to see which model reproduces the N 2 H’^ ring 
better. Since the radiative transfer calculation is out of the scope of the present work, we 
compare the column density of N 2 H+ in our full-network models with the estimated values 
in TW Hya. 


In the N 2 H'^ observation bv lQi et al.l (1201311. the 1-sigm a detection limit corresponds to 
the N 2 H’'' column density of 2 x 10^^ cm“^. iQi et al.l (120131 1 constructed disk models to fit 


their observational data. In the models that can reasonably fit the observational data, the 
N 2 H+ column density at the inner edge of the ring ranges from 4 x 10^^ to 2 x 10^® cm“^. 
They also found that the column density contrast at the inner edge of the N 2 H+ ring is at 
least 20, and could be larger. 


We also calculated the optical depth of NoH+ {J = 4 — 3) line for a slab of H 2 and N 2 H+ 
gas using the Radex code (Ivan der Tak et al. 1120071) . Assuming the H 2 density of 10® cm“® 
and line width of 0.15 km s“^, the N 2 H’'' column density of 10^^ cm“^ corresponds to the 
optical depth of r = 0.60 and 0.49 for the gas temperature of 17 K and 30 K, respectively. 
In order to reproduce N 2 H+ ring, the c olumn density o f N 2 H+ should be at least higher than 


10 


12 


cm 


-2 


, which is consistent with iQi et al.l (1201311 . 


Radial distributions of column densities of CO, HCO’*', N 2 and N 2 H+ in the dark cloud 
dust model are plotted in Figure [5])a) and (c). The horizontal axis, R, is linear in (a), while it 
is logarithmic to highlight the inner radius in (c). The distribution of N 2 H’^ column density 
is rather flat. At f = 1 x 10^ yr, the N 2 H+ abundance has a peak in the midplane at i? ~ 200 
AU, which is outside the CO snow line. But N 2 H''' column density does not sharply drop 
inwards at the CO snow line (~ 125 AU), because N 2 H’'' is abundant in the disk surface even 
in the inner radius. At 9.3 x 10^ yr, the N 2 H+ is abundant only in the disk surface, and the 
N 2 H+ column density shows a peak at ~ 30 AU, which is inside the CO snow line. 


The column densities of molecules in the millimeter grain model are shown in Figure [5] 
(b) and (d). At f = 1 x 10^ yr, the N 2 H+ column density has a peak at the radius of ~ 40 
AU, and signihcantly decreases inwards at ~ 20 AU, which corresponds to the CO snow line 
in the disk model with millimeter grains. This model is thus in better agreement with the 
observation of TW Hya than the dark cloud dust model. At f = 9.3 x 10^ yr, the radius of 
N 2 H+ column density peak is shifted inwards due to the sink effect on CO and N 2 . 
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4. Analytical Solution for Molecular Ion Abundances 

In the molecular layer, the major ions are HCO’*', Hjj' and N 2 H+. In this section we 
derive analytical formulas for their abundances. The analysis helps us to better understand 
the spatial distributions of these species (Figure [2] and |3]). We will show that we can calculate 
HCO'*' and N 2 H+ abundances in a disk model, if the distributions of gas density, tempera¬ 
ture, ionization rate, and abundances of CO and N 2 are given. The formulas are useful in 
constraining the ionization rate from the observation of CO, HCO"*" and N 2 H’''. 


4.1. Electron 


Electron abundance, i.e. the ionization degree, is determined by the balance between 
ionization and recombination. Let us hrst assume that HCO+ is the dominant ion in the 
disk, for simplicity. A sequence of reactions starts with ionization of H 2 via cosmic-ray. 
X-ray, or decay of radioactive nuclei; 


H 2 ^ H+ + e 


H:: 


H 2 ^ H+ + H 
H+ + CO ^ HCO+ + H; 
HCO+ + e^H + CO. 


(7) 

( 8 ) 
( 9 ) 

( 10 ) 


We define the ionization degree a:(e) as n(e)/nH rather than n(e)/n(H 2 ); while the latter is 
the usual definition, molecular abundances are relative to hydrogen nuclei in our numerical 
calculation. Considering the neutrality, the ionization degree in molecular gas is given by 


x{e) = 


n(e) n(e) 


c 


flu 


2n(H2 


2A[X0p' 


( 11 ) 


H 


The ionization rate C is a sum of the rates by cosmic ray. X-ray (Figured]), and decay of 
radioactive nuclei. The rate coefficients of relevant reactions are listed in Table [T]in Appendix 

A. 


Now, it should be noted that the dominant ion varies within the disk. Although the rate 
coefficients of the dissociative recombination of molecular ions are mostly of the same order 
(~ 10“^ cm^ s“^), the values vary slightly among molecular ions. We thus adopt an iteration 
to evaluate the ionization degree. Initially, we assume that the reaction ffTOl) dominates in 
the recombination of electron in the gas phase, and calculate the ionization degree, which is 
used to evaluate the abundances of molecular ions (see the following subsections). Then we 
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re-calculate the equation flTT]) by replacing /ijjg] with the average of the recombination rate 
coefficients of molecular ions weighted by their abundances. 

In the midplane at i? < a few 10 AU, the gas density is so high that the grain surface 
recombination becomes more effective than the gas-phase recombination; 

HCO+ + G(-) ^H + CO + G(0), (12) 

where G(-) and G(0) represent a negatively-charged grain and a neutral grain, respectively. 
Then the ionization balance is described as 


Cn(H 2 ) = A[jQp(HGO+)n(e) fcG^^(G-)r^(HGO+) = AjXQp^(e) kGn{G-)n{e), (13) 


wher e fcc is fhe rate coefficient of the grain surface recombination of HGO^ (see e.g. lUmebavashi 
19831) . Assuming that most grains are negatively charged, which is valid when uh/C ^ 10^° 
cm“^ s (lUmebavashi Ill983l) . the electron abundances is calculated by solving this quadratic 
function. It should be noted that the equation flT^ is equivalent to the equation (ITT]) , when 
the gas-phase recombination is more effective than the grain-surface recombination. Thus 
we use the equation flT3|l rather than ffTTj) in the calculation of ionization degree in the whole 
disk, and is replaced by the weighted mean of the recombination rate coefficients in the 
iteration. 

In deriving the equations (ITT]) and da, we have assumed that the molecular ions are 
the dominant charge carrier. In the disk surface however, atomic ions such as G^ and S’*" 
are produced via photoionization and dominate over molecular ions (§3.2). In the transition 
layer from such an atomic-ion dominated (AID) layer to the molecular layer, many reactions 
contribute to the ionization balance, and the major reactions vary over the transition layer. It 
is thus difficult to derive analytical formula of electron abundances there. In the present work, 
we compare the electron abundance obtained by the equation flT^ with the calculation of 
full reaction network at each position in the disk. We adopt the latter, when it is twice larger 
than the former. As we have seen in the previous section, the position of the transition region 
depends sensitively on the grain-size distribution in the disk model, as well as on the assumed 
UV flux from the central star and the interstellar radiation helds. Readers are advised to 
use the PDR (photon-dominated region) codes to calculate the electron abundance in the 
AID layer and transition region for a specihc disk model. The PDR codes or their analogues 
are commonly used to calculate the vertical distributions of gas density and tempe ratures 

(e.g. 


Kamp &: Dullemond Il2004 iNomura &: Millar 1 120051: iGorti &: Hollenbach 1120081) . 
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4.2. H+ 

Hjj" is produced by the reaction ([ 8 ]). The major destruction paths of Hjj" are recombi¬ 
nation and proton transfer to CO and N 2 . In the midplane at i? < a few 10 AU, the grain 


surface recombination also becomes effective; 

H^ + e^Hs + H or H + H + H (14) 

H+ + G(-) ^H 2 + H + G(0) or H + H + H + G(0) (15) 

H+ + GO ^ HGO+ + H 2 , (16) 

H 3 + + N 2 ^ N 2 H+ + H 2 , (17) 


by 


Gonsidering the balance between the formation and destruction, abundance is given 


x(H3 +) 


1 _ C/^H _ 

2 /|np(e) + A[T 5 p(G-) + /:(i 6 p(GO) + /fi7p(N2) ' 


(18) 


4.3. N 2 H+ 


N 2 H''' is formed by the reaction (IT7)l . and is destroyed by the proton transfer to GO 
(reaction [T]) and recombination in the gas phase and on negatively charged grains: 

N 2 H+ + e^NH + N or N 2 + H (19) 

N 2 H+ + G(-) ^NH + N + G(0) or N 2 + H + G(0). (20) 


Then its abundance is given by 

x(N2H+) 


_ A[T7p(H+)a:(N2) _ 

A{X9p(e) + A[ip(GO) + ' 


( 21 ) 


4.4. HCO+ 

HGO’*' is formed by the reaction ([T]) and fllOp . and destroyed by the recombination in 
the gas phase flTUD and on grain surfaces (IT^ . Once the dust temperature exceeds ~ 100 K, 
molecules with higher proton affinity than GO, such as NH 3 , are desorbed to the gas phase 
to destroy HGO"*". In this section, we neglect such a high temperature region, which is rather 
limited in our disk model. It is straightforward to derive 


ri(HCO+) top(Ht) + %T=(N2H+) 
n(CO) + l|i2fi(G-) 


( 22 ) 
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At the disk surface region, where the grain surface recombination is not effective and 
dominates over N 2 H+, the equation is modihed to 


n(HCO+) _ f^jCO) 


It shows that HCO’*' is less abundant than if the abundance ratio of CO to electron is 
lower than ^ = 1.5 x 10^{T/300K)-^-^^ (cf. §3.1). 


4.5. Comparison with numerical results 

We now check how the analytical formulas compare with the full network calculation. We 
adopt the physical parameters (density, temperature, and ionization rate) and abundances 
of CO and N 2 from the full network model at f = 1 x 10^ yr and calculate the abundances 
of electron, H^, N 2 ll^, and HCO’*' using the analytical formulas at each position in the disk 
models with dark cloud dust and millimeter grains. As described in §4.1., our analytical 
solution applies to a layer in which the major ions are molecular ions. When the electron 
abundance obtained in the full network calculation is twice larger than given by the analytical 
formula (i.e. AID layer), we adopt the electron abundance from the former. It should be 
noted however, that the analytical formulas for Hjj’, N 2 ll^ and HCO"*" are appropriate even 
in the AID layer, if the electron abundance is adopted from the full network. 

Figure |6] shows the 2D distributions of electron, Hjj', N 2 II+ and HCO’*' abundances 
calculated using the analytical formulas in the disk models with dark cloud dust (left panel) 
and millimeter grains (right panel). The distributions of molecular ions are to be compared 
with those in the left columns of Figure |2] and El The dashed lines indicate the height {Z) 
above which we adopt the electron abundance of the full network calculation. For a more 
quantitative comparison. Figure [7] shows the molecular distribution in the Z-direction at 
R = 53.4 AU in the disk models with dark cloud dust and millimeter grains. The solid lines 
depict the abundances calculated by the full reaction network, while the dotted lines depict 
the analytical solution. 

We can see that the analytical formulas are in reasonable agreement with the results of 
the full network calculation in the both models. At high Z (e.g. Z > 0.4 at i? ~ 30 AU), 
is overestimated; while the analytical formula fllSp assumes that the hydrogen is all in 
H 2 , it is photodissociated in the full network model in such un-shielded low-density regions. 
In molecular layers at lower Z, on the other hand, the analytic formulas tend to slightly 
overestimate the molecular ion abundances, partly because the analytical formulas neglect 
the neutral species other than CO and N 2 . In the full network model, there are minor neutral 
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species, such as OH, which have a larger proton affinity than N 2 and CO. In the upper layers 
of the millimeter grain model, the analytical formula underestimates the HCO’*' abundance; 
another formation path, CO’*' + H 2 becomes effective in this region. 

Figure [8] shows the column densities of HCO’*' and N 2 H~*" obtained by the full network 
calculation in the millimeter grain model at f = 1 x 10® yr (solid lines) and the analytical 
formulas (dotted lines). The column densities of the analytical model agree with the full 
network model within a factor of 2. Although the analytical formula underestimates the 
HCO'*' abundance at the disk surface, the surface region does not contribute much to the 
HCO'*' column density. 


In summary, we have demonstrated that the analytical formulas of molecular ion abun¬ 
dances agree well with the full network results. The analytical formulas give reasonable 
spatial distributions of molecular ions, taking into account the spatial variation of gas den¬ 
sity and ionization rate. They can thus be used to compare a disk model with molecular line 
observations without performing the full chemical network calculations. The input parame¬ 
ters of the formulas are density, temperature, ionization rate, and abundances of CO and N 2 . 
The distributions of density and temperature should be prepared for a specific object. Then 
the ionization rate can be calculated by X-ray radiation transfer and by assuming a c osmic- 


ray penetration depth and abundances of radioactive nuclei (e.g. ICleeves et ah II2014I) . The 


electron abundance in the disk surface is also needed, and can be obtained by the PDR 
calculation. Alternative option is to simply assume a height from the midplane below which 
molecular ions are more abundant than atomic ions. Finally we need spatial distributions of 
CO and N 2 . While the full network calculations show that they are subject to the sink effect 
and thus could decrease with time even in the region warmer than their sublimation tem¬ 
perature, the simplest assumption would be the equilibrium abundance between adsorption 
and desorption, as we will see in the next section. 


5. Results: No Sink Model 

In our full network calculations, N 2 H’'' ring is reproduced in a millimeter grain model 
with a significant CO depletion due to the sink effect. Efficiency of the sink effect, however, 
depends on various parameters. It is therefore useful to calculate the distribution of N 2 H+ 
in a model without the sink effect. As described in §2.3, here we assume that the total (gas 
and ice) abundances of CO and N 2 are constant, and that their gas/ice ratios are determined 
by the equilibrium between the adsorption onto and desorption from the grain surfaces (eq. 
m and [5]) . 
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A combination of the analytical formulas of molecular ions and equilibrium abundances 
of CO and N 2 make it very easy to investigate the dependence of N 2 H+ abundance on various 
disk parameters. Here we demonstrate this merit of the analytical formulas by investigating 
the dependence of N 2 H+ abundance on desorption energies of CO and N 2 , and on ionization 
rate in the disk model with millimeter grains. We could apply the no-sink model to the disk 
model with dark cloud dust, as well, but the relatively high abundance of N 2 H'*' in the disk 
surface makes the radial distribution of N 2 H''' column density flatter than observed (§3.1). 


5.1. N 2 H+ in Disk Models Without Sink Effect 

Figure [9] (a — c) shows the distributions of CO, N 2 , and N 2 H”^ in the no-sink model in 
the disk model with millimeter grains. The desorption energies of CO and N 2 are assumed 
to be 1150 K and 1000 K, respectively. The dashed line depicts the height above which 
the electron abundance is adopted from the full-network model. Figure M {d — e) shows 
the vertical distributions of molecules at the radius of 40 AU and 93 AU. At i? = 40 
AU, the temperature is close to the sublimation temperatures of CO and N 2 even in the 
midplane. The gas-phase abundances of CO and N 2 are thus mostly determined by the 
thermal desorption. In the midplane at i? = 93 AU, on the other hand, the non-thermal 
desorption by cosmic-ray dominates over the thermal desorption. The abundance of N 2 H”^ is 
determined by the equation fl?ID . As expected, N 2 H+ abundance is high in layers where N 2 is 
more abundant than CO. We note however, that N 2 H+ abundance is determined not only by 
the abundance ratio of N 2 /CO. For example, in Figured (d — e), the N 2 H’^ abundance varies 
among the three positions where the abundances of CO and N 2 are equal. The absolute 
values of CO and N 2 abundances matter, since they control the abundance of Hg (eq. ITS]) , 
from which N 2 H’'‘ is formed. Note that N 2 H''’ abundance in the very surface region of the 
disk (e.g. Z > 0.3 at i? ~ 50 AU) could be overestimated; its mother molecule, N 2 , is 
photodissociated in the disk surface in the full network model. 

The radial distributions of N 2 H~'“ column density is shown in Figure O (/). The dot¬ 
ted line depicts the CO column density multiplied by a factor of 10“^. In order to avoid 
the photodissociation region at the disk surface, N 2 H+ column density is calculated by the 
integration at \Z\ < 0.3 (R > 50 AU), \Z\ < 0.15 (10 < i? < 50 AU), and \Z\ < 0.1 
{R < 10 AU), while the CO column density is obtained by the integration along the whole 
disk height. We can see that N 2 H''' column density indicates a ring structure, as in the 
full-network model. The spatial distribution of N 2 H’^ is however, different from that in the 
full-network model. Firstly, the N 2 H+ ring is sharp and exists right outside the CO snow line 
in the no-sink model, simply reflecting the lower desorption energy of N 2 than that of CO, 
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while it is broader in the full-network model, especially in the late stage {t = 9.3 x 10^ yr). 
The peak N 2 H”'' column density at i? ~ 50 AU and t > 10^ yr is lower in the full network 
model, because N 2 is depleted via the sink effect in the midplane. Secondly, N 2 H’'' is conhned 
to a thin layer colder than the CO sublimation temperature in the no-sink model, while it 
is abundant also in the upper warmer layers in the full-network model. The two models can 
thus be distinguished by constraining the vertical distributions of N 2 H''' from observations. 


I n the case of T W Hya, N 2 H’^ lines of J = 3 — 2 and J = 4 — 3 are observed flQi et la. 


2 OI 3 I: lOi et al.l 120131) . Assuming that the lines are optically thin under LTE conditions, the 


excitation temperature is derived to be 35 ±10 K. The relatively high excitation temperature 
of N 2 H'*~, togethe r with the low CO abundance indicated by C^®0 and HD observations 
flFavre et al.l 120131) . might be better explained by the full-network model. We postpone the 
discussion on TW Hya to a future work, in which we will apply our network model and 
analytical formulas to a disk model specified for TW Hya, and simulate molecular emission 
with non-LTE radiative transfer calculations for more quantitative comparisons with the 
observational data. 


5.2. Dependence on desorption rates 


So far, we set the desorption energies of CO and N 2 to be 1150 K and 1000 K, respec¬ 
tively. It is well known however, that the desorption energies of molecules depend on the 
chemical compositions and physical st ructure (e.g. cryst a l or a morphous) of the ice man¬ 
tle. While the values we adopted from iGarrod fc Herbst (120061) are desorption energies on 
water-domina ted ice surfaces deri ved from the temperature programmed desorption (TPD) 
experiments (jCollings et ahl 120041) . the desorption en ergies of pure ices of CO and N 2 are 
Edes(CO)= 855 ± 25 K and Edes(N 2 ) = 790 ± 25 K f Oberg et ah 2005). In §3.2, we have 
shown that in the full-network model the N 2 H+ distribution does not sensitively depend on 
the desorption energies of CO and N 2 , because the layer of CO depletion, where N 2 H+ is 
abundant, is determined by the sink effect rather than thermal desorption. If the CO sink 
is not effective, on the other hand, the abundance of N 2 H''' could be more sensitive to their 
desorption energies. 

We also note that the equation (jl]) assumes the first-order desorption; i.e. the desorption 
rate is proportional to the abundance of the species in the ice mantle, ncoice- In the following, 
we call this model as ”2 phase”, since this equation is usually used in the 2 -phase gas-grain 
chemical models, which do not discriminate the layers of ice mantle. In reality, the migration 
(and thus desorption) of molec ules deeply embedded in the ice mantle could be inefficient, at 


least at low temperatures (e.g. ICollings et al.ll2004l) . We therefore consider another model. 
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i.e. ” 3-phase” model, in which only the surface monolayer is subject to desorption: 


S'Tra^Uth’^cogas^^dust = min[47ra^ndustA^site,?^coice] j i^exp 


kT ) 


+ VTcKCpe exp 


E2^ 

kTyrtPiM 


(24) 

where N^ite = 1-5 x 10^® cm“^ is the number density of adsorbing site on a grain surface. 
The desorption is zero-th order, as long as the ice is abundant enough to occupy more than 
a monolayer in the ice mantle, i.e. ricoice > dvra^Udust-^site- Then the gaseous abundance is 
independent of the total abundance. In the 2-phase model, on the other hand, the gaseous 
CO (N 2 ) abundance is proportional to the assumed total abundance of CO (N 2 ). 


Before showing the N 2 H’^ abundance in the disk model with various desorption rates, 
it is instructive to apply the analytical formulas to a simpler model, where the gas density 
and ionization rate are constant, 10® cm“® and 5 x 10“^^ s“l, respectively. Figure fTOl shows 
the gaseous abundances of CO, N 2 and N 2 H''’ as a function of temperature. The gas and 
dust temperatures are set to be equal. The total (gas and ice) abundances of CO and N 2 
are 1 x 10“^ and 4.5 x 10“®, respectively. The desorption energies are set to be F^des(CO) = 
F^des (^ 2 )= 855 K in Figure fTOfah and F^des(CO)= 855 K and F^des(N 2 )= 790 K in Figure 
[TOlf bh The desorption rate of 2-phase model (eq. H]) is assume for the red lines, while 
the 3-phase model is assume for the blue lines. The abundances change drastically at the 
sublimation temperature ~ 20 K. At T < 15 K, the gaseous abundances of CO and N 2 
slightly increase with decreasing temperature; desorption rate is kept constant due to the 
non-thermal desorption, while the adsorption rate is proportional to By comparing the 
red and blue lines, we can see that the sublimation temperature and gaseous abundances of 
CO and N 2 at low temperatures are signihcantly different between the 2-phase and 3-phase 
models; the 3-phase model gives much smaller desorption rate than the 2-phase model, when 
the ice mantle is thick. In the 3-phase model with the same E^es for CO and N 2 (blue 
lines in panel a), the gaseous abundance of CO and N 2 are the same below the sublimation 
temperature. It is interesting that the N 2 H+ abundance shows a sharp peak slightly below 
the sublimation temperature. It clearly shows that the N 2 II+ abundance depends not only 
on the abundance ratio of CO/N 2 , but also on the absolute values of their abundances. 
Specihcally, the N 2 H''' abundance reaches the maximum value when the abundance ratio of 
CO to electron is ~ which is about 3.3 x 10® at T = 17 K (see Appendix B). In 

the models of Figure [TUI the electron abundance is about several xl0“®° and varies slightly 
with temperature and desorption model. At warm temperatures where the abundance ratio 
of CO to electron is higher than this critical value, HCO’*' is the dominant ion, while 
dominates at lower temperatures. In the case of 2-phase model with the same for CO 
and N 2 (red lines), on the other hand, N 2 II®' abundance does not show such a sharp peak 
around the CO sublimation temperature, because the abundance ratio of CO to electron is 
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always higher than the critical value mentioned above. In the case of E^es (N 2 ) < E^es (CO) 
(Figure [To] b), there is a narrow temperature range in which N2 is relatively abundant but 
CO is not, so that N 2 H+ abundance has the maximum value, even in the 2-phase model. 

Figure [TT] shows the radial distributions of N 2 H+ (solid lines) and CO (dotted lines) 
column densities in the disk model with millimeter grains as in Figure 9 (/), but for three 
sets of desorption energies of CO and N 2 : i?des(CO)=1150 K and i?des(CO)=1000 K in 
panel (a), i?des(CO)=855 K and i?des(CO)=790 K in panel (b), and i7des(CO)=855 K and 
F'des(CO)=855 K in panel (c). The column density of hydrogen nuclei is multiplied by a factor 
of 10“^^, and is shown by the green line. The black and red lines depict the 2-phase and 
3-phase models, respectively. The total (gas and ice) abundance of N 2 is 4.5 x 10“® relative 
to hydrogen for the solid lines. In order to investigate the dependence of N 2 II+ abundance 
on the total N 2 abundance, it is set to be 10 times higher for the dashed lines. We can see 
that all the models show a ring-like structure of N 2 II’*’. The N 2 H^ column density and its 
radial gradient however, depend significantly on the desorption energies and the model of 
desorption. The CO snow line and the inner edge of N 2 II’*’ ring are at smaller radii in the 
model with higher i7des(CO). They are also at smaller radii in the 3-phase model than in the 

2- phase model, because of the lower desorption rate in the 3-phase model. While the N 2 H+ 
column density signihcantly depends on the total N 2 abundance in the 2-phase model, the 
dependence is weak in the 3-phase model. Dependence of the N 2 II”*' column density on the 
desorption energies of CO and N 2 are complex, as expected from Figure [TOl In the models 
with i7des(CO) > i^des (N 2 ), there is a region where N 2 is thermally desorbed but CO is not, 
which results in a peak of N 2 H’'‘ column density with a width of a few 10 AU. In the models 
with Edes{CO) = Edes (N 2 )= 855 K, the N 2 H+ column density shows a sharp peak in the 

3- phase model, but not in the 2-phase model, as expected from Figure [TOl 


5.3. Dependence on Ionization rate 


Al though we have assumed that the attenuation length of the cosmic ray ionization is 96 
g cm~^ fjUmebavashi fc Nakano Ill98ll) . the stellar winds and/or magnetic fields could hamper 


the penetration of co s mic rays to the disk. T hen X-ray would be the major ionization source 
(jCleeves et al. II2013I: iGlassgold et al.lll997n . The dotted lines in Figure [3] depict the height 
from the midplane below which the X-ray ionization rate is lower than 10“^^ s“^ and 10“^® 
s“^, respectively. The layer with abundant N 2 H''' are around or below these lines in both 
the full network model (i.e. with CO sink) and the no-sink model. Here we investigate how 
the N 2 H’^ column density changes, if the cosmic ray does not reach the disk. 


Figure [T2I shows the radial distribution of N 2 H+ column density without the cosmic- 
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ray ionization. X-ray is the only ionization sonrce for the dashed lines, while the decay of 
radioactive nuclei is considered for the solid lines. It should be noted that cosmic ray also 
causes non-thermal desorption in our models; the non-thermal desorption by the cosmic-ray 
is neglected for the solid lines and dashed lines. For a comparison, the dotted lines depict the 
model in which the non-thermal desorption is included but ionization source is X-ray and 
radioactive nuclei; although such a model is not self-consistent, we can see the importance 
of the non-thermal desorption by comparing the solid lines with dotted lines. The blue and 
green lines depict the model in which CO and N 2 abundances in the gas phase are given by 
the 3-phase model; labels in the hgure depict the assumed desorption energies of CO and 
N 2 . For the red lines, we refer to the CO and N 2 abundances in our full network model at 
1 X 10^ yr, i.e. the model with sink. Note that there is thus no red dotted line. Compared 
with the models with cosmic ray ionization (Figure [H] and [TT]), the N 2 H’'' column density 
is reduced signihcantly. The reduction factor, the ratio of the peak column density in the 
model with cosmic ray to that in the model without cosmic-ray, is the highest (~ 70) in the 
no-sink model with low CO and N 2 desorption energies, while it is the lowest (factor of ~ 20) 
for the model with sink. In the latter model, N 2 H’'' is relatively abundant even in the layer 
with the X-ray ionization rate > 10“^® s“^ (i.e. Z/i? > 0.2). In the no-sink models, on the 
other hand, N 2 II''' is depleted at 0.2 < Z < 0.3. In the model with lower desorption energy 
of CO, the CO freeze-out region, where the N 2 ll’^ is abundant, is conhned to the layer closer 
to the midplane with low X-ray ionization rate. 


6. Summary 


In this work we calculated the molecular abundances in disk models to invest igate the 


origin of N 2 H’'‘ ring recently found at the CO snow line in the disk of TW Hya (jOi et ah 


20131 ). We adopt two disk models with different dust properties; dust grains are similar 


to dark cloud dust in one model, while they have grown to millimeter size in the other 
model. We hrst calculated a full network of gas-grain chemical reactions. Our hndings in 
the full-network model are as follows. 


• In the model with dark cloud dust, N 2 II"'' is abundant in the disk surface. Although 
N 2 H''' abundance near the midplane has a local maximum outside the radius of CO 
snow line, N 2 H’'‘ column density is rather high even inside the CO snow line, because 
N 2 H''' in the disk surface contributes signihcantly to the column density. 

• In the model with millimeter grains, the column density of N 2 II+ shows a local peak 
around the CO snow line. N 2 II+ is abundant in the warm intermediate layer where 
CO is depleted via the sink effect, while N 2 II'*’ is destroyed by the proton transfer to 
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CO inside the CO snow line. Penetration of UV radiation (i.e. lower extinction than 
the dark clond dust model) plays two important roles. Firstly, UV radiation makes 
the atomic ions dominant in the disk surface, so that N 2 H+ is not abundant there. In 
the intermediate layer where CO is depleted via the sink effect, the conversion of N 2 to 
NH 3 (i.e. the sink effect on N 2 ) is prevented by the photodissociation of NH, so that 
N 2 H''' is abundant. 

• In the model with millimeter grains, the distributions of N 2 H~*' and its column density 
do not change signihcantly when the desorption energies of CO and N 2 are varied. 
The N 2 H'^ column density shows a local peak at the radius of CO snow line, even if 
the desorption energies of CO and N 2 are equal, because in the region with abundant 
N 2 H''', CO is depleted via the sink effect rather than the adsorption of CO itself onto 
grain surfaces. 

• By analyzing the results of the full-network model, we derived analytical formulas of 
electron, H 3 , N 2 H+ and HCO+ abundances as functions of gas density, temperature, 
ionization rate, and abundances of CO and N 2 . The analytical formulas would be 
useful for radio observers to derive the abundances of these molecular ions from the 
observational data using a reasonable physical model for a well-observed disk such as 
TW Hya. 

While the sink effect on CO plays an important role in determining the N 2 H’'' abundance 
in the full-network model, the efficiency of the sink depends on various parameters such as 
turbulence in the disk and the rates of chemical conversion of CO to less volatile species. We 
thus constructed the no-sink model, in which the total (gas and ice) abundances of CO and 
N 2 are set constant, and their gas-phase abundances are determined by the balance between 
the adsorption and desorption. The abundances of molecular ions are calculated using the 
analytical formulas. The results of the no-sink models are as follows. 

• The column density of N 2 H~'' shows a ring-like structure in the no-sink model with 
millimeter dust grains. Since the abundance of N 2 H+ is given by a non-linear function 
of CO and N 2 abundances, it can reach a moderate value even in the cold N 2 freeze-out 
region, depending on the abundances of CO and N 2 . 

• Even if the desorption energies of CO and N 2 are the same, N 2 H'^ abundance peaks at 

the temperature slightly below the CO (and N 2 ) sublimation temperature, where the 
abundance ratio of CO to electron is ~ ~ 10^- 

• Although the N 2 H+ ring is produced both in the full network model and in the no-sink 
model, the detailed distributions of N 2 H+ are different between the two models. In 
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the no-sink model, N 2 H~'“ abundant layer is confined to a layer colder than the subli¬ 
mation temperature of CO, while in the full-network model (i.e. with CO sink), N 2 H+ 
abundant layer extends to warmer layers. These models can thus be discriminated in 
observations if we can determine the vertical distribution of N 2 H+ in disks, specifically, 
starting with constraining the excitation temperature of N 2 H+ emission. 

• The column density of N 2 H+ in the no-sink model sensitively depends on the desorption 
rates of CO and N 2 . 

• If the cosmic ray does not reach the disk, the N 2 H'^ column density is reduced by a 
factor of 20 in the model with CO sink, and by a factor of 70 in the no-sink model. 
The reduction is more significant in the model with lower desorption energy of CO, 
because N 2 H'^ (i.e. CO freeze-out) is confined to a layer closer to the midplane, where 
X-ray ionization rate is lower. 
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A. Reaction Rate Coefficients 


We list the rate coefficients of the reactions relevant to the analytical formulas of molec¬ 
ular ions in Table [H In the full network model, the spin state (ortho and para) of is 
discriminated, and the reaction ra tes actually vary w i th the spin state. For those reactions, 
we use the rate coefficients from Gar rod &: Herbst fj2006 ') in the analytical formulas, for 
simplicity. 


B. Dependence of N 2 H+ Abundance on Temperature 

Figure dO] shows that N 2 H+ abundance reaches the maximum value at T ~ 17 K in 
the 3-phase model, even if we assume the same desorption energy for CO and N 2 . Here we 
analyze the model to show that N 2 H+ abundance is maximized when the abundance ratio 
of CO to electron is ~ 10^. 
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At the density of nn = 10® cm“®, gas-phase recombination is more effective than the 
grain surface recombination. Then the analytic formula of the N 2 H+ abundance (eq. [?I1) 
becomes 


^(N^H+) = , fiy(H3 )^(N2) 


<(I 2 p(e) + *lip(CO) ' 

When the abundance ratio of CO to electron is higher than which is about 3.3 x 

10® at T = 17 K, the denominator is dominated by the second term. Substituting the 
analytical formula of Hjj" abundance (eq. [T5]) to (IBip . it is straightforward to show that 
N 2 H''‘ abundance is proportional to x(N 2 )/a;^(CO), and thus increases with decreasing CO 
abundance at T ~ 17 — 20 K. If the hrst term dominates in the denominator, on the other 
hand, N 2 H’^ abundance is proportional to N 2 abundance. The N 2 H’'' abundance decreases 
with decreasing temperature at T ~ 15 — 17 K. 
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Table 1: Reaction rate coefficients relevant to the analytical formulas®" 


reaction 

a 

/3 

H+ + CO ^ HCO+ + Ha 

1.61(-9)^ 

0.0 

H+ + Na ^ NaH+ + Ha 

1.7(-9) 

0.0 

H+ + e ^ Ha + H 

2.59(-8) 

-0.5 

H+ + e^H + H + H 

4.61(-8) 

-0.5 

HCO++e^H + CO 

2.40(-7) 

-0.69 

NaH+ + CO ^ HCO+ + Na 

8.8(-10) 

0.0 

NaH+ + e ^ NH + H 

1.30(-8) 

-0.84 

NaH+ + e ^ Na + H 

2.47(-7) 

-0.84 


“Rate coefficients are given in the format k = a x (T/300.0)^ cm^ s 
'^A{B) stands for A x 10^ 
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Table 2: Elemental abundance and initial abundances of assorted molecules 


Element 

Abundance'^ 

Element 

Abundance 

H 

1.0 

He 

9.75(-2) 

N 

2.47(-5) 

0 

1.80(-4) 

C 

7.86(-5) 

S 

9.14(-8) 

Si 

9.74(-9) 

Ee 

2.74(-9) 

Ee 

2.74(-9) 

Na 

2.25(-9) 

Mg 

1.09(-8) 



Species 

Abundance 

Species 

Abundance 

H 2 O 

1.15(-4) 

CO 

3.57(-5) 

CO 2 

3.52(-6) 

CH 4 

1.50(-5) 

H 2 CO 

1.20(-5) 

CH 3 OH 

6.50(-6) 

N 2 

4.47(-6) 

NHs 

1.44(-5) 


A{—B) means T x 10 ^. 
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Fig. 1.— Distribution of the density of hydrogen nuclei nH(top), gas temperature, dust 
temperature, and ionization rate by X-rays (x (bottom) in the disk model with dark cloud 
dust (left panels) and millimeter grains (right panels). 
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Fig. 2.— Distributions of CO, HCO’*', N 2 , N 2 H+, and electron in the gas phase at 
t = 1 X 10® yr (left panels) and 9.3 x 10® yr (right panels) in the model with dark cloud 
dust. Dashed lines in the panels of CO and N 2 depict the position where the gas-phase and 
ice-mantle abundances become equal in the adsorption-desorption equilibrium (eq. [6]). The 
dotted lines depict the position where the X-ray ionization rate is equal to the cosmic-ray 
ionization rate (5 x 10“^^ s“^) and to the ionization rate by decay of radioactive nuclei 
(1 X 10-^8 
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Fig. 2.— cont. 
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Fig. 3.— Distributions of CO, HCO’*', N 2 , N 2 H''’, Hg and electron as in Figure 2, but for 
the model with millimeter grains. 
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Fig. 3.— cont. 
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Fig. 4.— Distributions of CO, HCO’*', N 2 , and N 2 H’'' in the model with millimeter grains as 
in Figured The desorption energies of CO and N 2 are set to be equal (1000 K). 
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Fig. 5. — Radial distributions of column densities of CO, HCO’*', N 2 and N 2 H''’ in the gas 
phase at t = 1 X 10^ yr (solid lines) and 9.3 x 10^ yr (dotted lines) in the model with dark 
cloud dust (a and c) and millimeter grains (b and d). The column densities of CO and 
N 2 are multiplied by a factor of 10“^ to £t in the hgure. The horizontal axis, R, is linear 
in (a) and (b), while it is logarithmic to highlight the inner radius in (c) and (d). The 
dashed lines in the right panel depict the column densities of molecules in the model with 
Edes(CO)= Edes(N 2 ) = 1000 K at 9.3 x 10® yr. 
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Fig. 6.— 2D distributions of HCO’*', N 2 H~'', Hjj", and electron calculated using the analytical 
formulas in the disk models with dark cloud dust (left panel) and millimeter grains (right 
panels). CO and N 2 abundances are adopted from the full-network model at f = 1 x 10® yr. 
The dashed line depicts the layer above which the electron abundance is adopted from the 
full network model. 
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Fig. 7.— Molecular abundances as a function of height from the midplane {Z/R) aX R = 53.4 
AU in the disk models with dark cloud dust (left panel) and millimeter grains (right panel) 
at t = 1 X 10^ yr. The solid lines depict the results of the full-network model, while the 
dotted lines depict the abundances obtained by the analytical formulas. 
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Fig. 8.— Radial distributions of column densities of HCO’*' and N 2 H+ in the full-network 
model with millimeter grains at f = 10^ yr (solid lines). The column densities obtained by 
the analytical formulas are depicted by the dotted lines. 






39 




Fig. 9.— (a-c) Distributions of abundances of CO, N 2 and N 2 H^ in the no-sink model. The 
desorption energies are set to be ii^des(CO)= 1150 K and F^des(N 2 ) = 1000 K. The dashed 
lines depict the layer above which the electron abundance is adopted from the full-network 
model, (d-e) Vertical distributions of molecules at i? = 40 AU and 93 AU. (f) Column 
densities of N 2 H’^ and CO. The column density of CO is multiplied by a factor of 10“^ to £t 
in the figure. 
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Fig. 10.— Abundances of CO (dotted lines), N 2 (dashed lines) and N 2 H+ (solid lines) as 
a function of temperature. The total (gas and ice) abundances of CO and N 2 are assumed 
to be 1 X 10“^ and 4.5 x 10“®, respectively. The gas density is set to be 1 x 10^ cm“^ and 
ionization rate is 5 x 10“^^ s“^. As for the desorption rate, the 2-phase model (eq. H]) is 
assumed for the red lines, while the 3-phase model (eq. [2l]) is assumed for the blue lines. 
The desorption energies of CO and N 2 are 855 K in panel (a), while i^des (CO) is 855 K and 
h^des (N 2 ) is 790 K in panel (b). 



Fig. 11.— Column density of N 2 H’'' in the no-sink models; the total (gas and ice) N 2 
abundance is set to be 4.5 x 10“® for the solid lines, while it is 4.5 x 10“® for the dashed 
lines. Black and red lines represent the 2-phase and 3-phase models, respectively. The dotted 
lines depict the column density of CO, which is multiplied by a factor of 10“^ to £t in the 
figure. The green lines depict the column density of hydrogen nuclei multiplied by a factor 
of 10“^^. N 2 H’^ column density is calculated by the integration at Z < 0.3 {R < 50 AU), 
Z <0.15 (10 < i? < 50 AU), and Z < 0.1 (i? < 10 AU) to avoid the photodissociation region 
at the disk surface, while CO and hydrogen column densities are obtained by integrating 
over the whole disk height. 
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Fig. 12.— Column density of N 2 H'*' in the models without cosmic ray ionization. X-ray is the 
only ionization source for the dashed lines, while the decay of radioactive nuclei is considered 
for the solid lines and the dotted lines. The non-thermal desorption by the cosmic-ray is 
neglected for the solid lines, while it is included for the dotted lines. For the blue and green 
lines, gaseous CO and N 2 abundances are given by the 3-phase model with the desorption 
energies labeled in the hgure. For the red lines, CO and N 2 abundances are adopted from 
the full-network model at 1 x 10^ yr, in which the sink effect on CO is at work. 














